% close all
% clear all
%symbols
% i = patch = rows
% j = age = columns
% 

%initial conditions
% Nij=Nij_initial; %Initial number of fish. rows=patches; cols=size classes
% Bij=Nij.*Wij; %biomass in each patch and age class

Bij=Bij_initial;
Nij=Bij./Wij;

%% simulation
tolerance=1000; Kdiff=tolerance*2;
%set up empty data matrices
Bijy=NaN(numpatches,max_age,max_simulation_length); %for recording time series of age-structured biomass
Nijy=NaN(numpatches,max_age,max_simulation_length);
Yiy=NaN(numpatches,max_simulation_length); %for recording time series of spatial yield
Revenue_iy=NaN(numpatches,max_simulation_length); %for recording time series
Cost_iy=NaN(numpatches,max_simulation_length); %for recording time series
Profit_iy=NaN(numpatches,max_simulation_length); %for recording time series

x0_aveProfit_start=1e8; %1e8  THIS EVENTUALLY WORKED FOR EQUIL. IT SEEMS TO TAKE FOREVER FOR NPV. TRY BELOW:
x0_aveProfit_start=x0_aveProfit_start/5;
% x0_aveProfit=10000; %initial guess (in future time steps update with previous choice)
% x0_aveProfit=1e4; %TargetBiomass=1e6
% x0_aveProfit=1e6; %TargetBiomass=1e8
x0_aveProfit=x0_aveProfit_start; %TargetBiomass=1e10
% energy_locs=zeros(numpatches,1); %start with no energy facilities
y=0;
beta_i=beta_i_initial;
while y<max_simulation_length
% while Kdiff>tolerance
    y=y+1;
    
    Kold=sum(Bij(:)); 
    
    Bi_mature=sum(Bij(:,age_mature:end),2); %biomass of mature fish in each patch
    Li=Bi_mature; %number of larvae produced in each patch ~ fish biomass
    Si=(Li'*Dii)'; %disperse the larvae according to transition probabilities
    RNi=(Si.*alpha1)./(1+Si.*beta_i); %number of recruits in each patch following DD
    Bi_legal=sum(Bij(:,age_legal:end),2); %biomass of legal size (set by age) fish in each patch
%     %%fleet model%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Fsum>0
%     Bi_legal_fish=Bi_legal; 
%     Bi_legal_fish=Bi_legal.*(energy_locs==0); %fishery can't take from energy sites

    Bi_legal_fish=Bi_legal.*(MPA_locs==0)'; %fishery can't take from energy sites
    
%     if sum(Bi_legal_fish)>0
        %rearrange this:
%         theta=pr*((Kvirgin/numpatches)*theta_percent)*(1-exp(-delta))*(1/(delta))
%         so this: (theta*delta)/((1-exp(-delta))*pr)=((Kvirgin/numpatches)*theta_percent)
    if sum(Bi_legal_fish)>(theta*delta)/(pr*(1-exp(-delta)));    
        
        [aveProfitfzero,Fdiff,EXITFLAG]=fzero(@Fleet_v1,x0_aveProfit); %,optionsfzero);
        %if Fsum is too high combined total habitat too low, then neg profit will result. Set to zero profit
        if aveProfitfzero<0; aveProfitzero=0; EXITFLAG=1; 
%             disp('Fsum negative, set aveProfitzero=0 and EXITFLAG=1'); 
        end
        while EXITFLAG==-6 
            disp('EXITFLAG==-6')
%             break
            x0_aveProfit=x0_aveProfit/5; %or divide by 5? 10?
            if x0_aveProfit<x0_aveProfit_start/1e5
                disp('Fleet cannot solve and x0_aveProfit shrinking toooooo small........break')
                break
            end
            [aveProfitfzero,Fdiff,EXITFLAG]=fzero(@Fleet_v1,x0_aveProfit); %,optionsfzero);
             	%if Fsum is too high combined total habitat too low, then neg profit will result. Set to zero profit
                 if aveProfitfzero<0; aveProfitzero=0; EXITFLAG=1; disp('Fsum negative, set aveProfitzero=0 and EXITFLAG=1'); end
       
            if EXITFLAG<=1; disp(['Problem solved with lower x0_aveProfit_start: EXITFLAG=',num2str(EXITFLAG),]); end
        end
        x0_aveProfit=aveProfitfzero; %use for next time step
        EXITFLAG_store(y)=EXITFLAG;
    else
        aveProfitfzero=0; %no harvestable fish, so no profit
        EXITFLAG_store(y)=0;
    end
        aveProfitfzero;
        
        

 
    
%     Fi_vector = (Bi_legal_fish.*pr-aveProfitfzero-aveProfitfzero*delta-theta-delta*theta)./(aveProfitfzero+theta); %(INVERSE MODEL)
%     Fi=max(Fi_vector,0);
       
    %(EXP MODEL)
    warning off MATLAB:divideByZero
    % Q=(aveProfitfzero+theta)./(pr+Bi_legal_fish); %pr PLUS Bi_legal_fish.......WRONG, SEE FleetModel_exp_v2.nb Mathematica code
    Q=(aveProfitfzero+theta)./(pr*Bi_legal_fish); %pr MULTIPLIED BY Bi_legal_fish.......RIGHT
    
    lambert_z=-exp(-1./Q)./Q;
    warning on MATLAB:divideByZero
    %wapr will break if lambert_z is out of range (i.e., <-exp(-1)); so set
    %those values to the edge of the range = -exp(-1)
    lambert_z(lambert_z<-exp(-1))=-exp(-1);
    lambert_w=wapr(lambert_z);
    d_plus_Fi=1./Q+lambert_w';
    d_plus_Fi=max([zeros(length(d_plus_Fi),1) d_plus_Fi],[],2); %get ride of tiny negative results
    Fi_vector=d_plus_Fi-delta;
    Fi=max(Fi_vector,0);   

elseif Fsum<=0 %Note: Fsum neg correction (
    Fi=zeros(numpatches,1);
end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    Survivei_legal=exp(-(Fi+delta)); %survival rate of harvested fish in each patch ~ fishing + natural mortality. Same for all legal-to-fish ages (EXP MODEL)
%     Survivei_legal=1./(1+Fi+delta); %survival rate of harvested fish in each patch ~ fishing + natural mortality. Same for all legal-to-fish ages (INVERSE MODEL)
    Yiy(:,y)=(Bi_legal.*(1-Survivei_legal)).*(Fi./(Fi+delta)); %biomass yield in each patch 
    Revenue_iy(:,y)=Yiy(:,y).*pr; %revenue in each patch
%     Cost_iy(:,y)=theta.*(log(Bi_legal)-log(Bi_legal-Yiy(:,y)));%cost in each patch
    Cost_iy(:,y)=theta.*Fi;%cost in each patch
    Profit_iy(:,y)=Revenue_iy(:,y)-Cost_iy(:,y); %profit in each patch
    
    Surviveij_legal=repmat(Survivei_legal,1,length(age_legal:max_age)); %survival rate per patch, per age class
    Survivei_nonlegal=exp(-(Fi.*0+delta)); %survival rate of non-harvested fish ~ natural mortality (EXP MODEL)
%     Survivei_nonlegal=1./(1+Fi.*0+delta);%survival rate of non-harvested fish ~ natural mortality (INVERSE MODEL)
    Surviveij_nonlegal=repmat(Survivei_nonlegal,1,length(1:(age_legal-1))); %survival rate per patch, per age class
    Surviveij=[Surviveij_nonlegal,Surviveij_legal]; %stitch survival rates together for all age classes
%     Nij_maxAgebin=Nij(:,end).*Surviveij(:,end); %survival of adults in final age class (survive and stay in that class)
    Nij(:,2:end)=Nij(:,1:(end-1)).*Surviveij(:,1:(end-1)); %number of fish that survive to next age class 
%     Nij(:,end)=Nij(:,end)+Nij_maxAgebin; %add to final age class adults in that class that survived
    Nij(:,1)=RNi; %add recruit number to the first age class
    
    Nijy(:,:,y)=Nij;
    Bij=Nij.*Wij; %convert # to biomass
    Bijy(:,:,y)=Bij;
    K=sum(Bij(:));  
%     Kdiff=abs(K-Kold);
   
end
%get rid of NaNs at end of records
%Not needed if sim length fixed (e.g, 100 years)
Yiy=Yiy(:,1:y);
Bijy=Bijy(:,:,1:y);
Nijy=Nijy(:,:,1:y);
Revenue_iy=Revenue_iy(:,1:y);
Cost_iy=Cost_iy(:,1:y);
Profit_iy=Profit_iy(:,1:y);

%scaled payoff - TEST
% tmp1=sum(Profit_iy,1);
% tmp2=tmp1(end); 
% payoff_scaled=tmp2/FleetProfit_Fsum_opt

%compensation ratio
% a=origin
% b=from origin to max number of settlers
%CR=(slope a)/(slope b)
% slope a (i.e, at S=0)
slope_a=alpha1;
% slope b (i.e., at S=Svirgin=Smax)
slope_b=sum(RNi)/sum(Si);
%note: also works for slope_b=sum(RNi(n))/sum(Si(n));
CR=slope_a/slope_b;
K; %sysem wide biomass

%% figures
if hide0_show1_figures==1
figure
plot(1:y,(squeeze(sum(Bijy,2)))')
xlabel('Year')
ylabel('Biomass per patch')
set(gcf,'color',[1 1 1]);box off

figure
plot(1:y,Yiy')
xlabel('Year')
ylabel('Yield per patch')
set(gcf,'color',[1 1 1]);box off

figure
plot(1:y,sum(Yiy))
xlabel('Year')
ylabel('Total yield')
set(gcf,'color',[1 1 1]);box off

figure
plot(1:y,Profit_iy)
xlabel('Year')
ylabel('Profit per patch')
set(gcf,'color',[1 1 1]);box off

figure
plot(1:y,sum(Profit_iy,1))
xlabel('Year')
ylabel('Total profit')
set(gcf,'color',[1 1 1]);box off

% figure
% plot(1:y,(squeeze(sum(Bijy,1)))')
% xlabel('Year')
% ylabel('Biomass per age class')
% set(gcf,'color',[1 1 1]);box off
    
figure
plot(1:y,(squeeze(sum((sum(Bijy,1)))))')
xlabel('Year')
ylabel('Total biomass')
set(gcf,'color',[1 1 1]);box off

% figure
% plot(1:y,(squeeze(sum((sum(Nijy(:,age_mature:end,:),1)))))')
% xlabel('Year')
% ylabel('Total abundance of adults')
% set(gcf,'color',[1 1 1]); box off

end

%if calculating NPV, then first determine it by multiplying each year by
%the discount rate and summing the years. 
%Then, put that solution in the final year value of the output Profit_iy
%and Bijy matrixes. We do this because the code pulls the final (supposedly
%equilibrium) year values for the analysis. But now the final values are
%not equilibrium but the NPV!
if NPVon==1
    tmp1=Profit_iy.*discounted_value_iy;
    tmp2=sum(tmp1,2); %sum the years
    Profit_iy=repmat(tmp2,1, max_simulation_length);
    
    tmp1=Bijy.*discounted_value_ijy;
    tmp2=sum(tmp1,3); %sum the years
    tmp3=NaN(numpatches, max_age, max_simulation_length); %shell
    for tmp4=1:length(tmp3(1,1,:));
        tmp3(:,:,tmp4)=tmp2;
    end
    Bijy=tmp3;
end
    

    
    
    
    
